library(survey)
library(ggplot2)

## run "PCM_cleaning.R"

## Create weighted survey design object
milconf_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

#### CONFIDENCE: Q10A BY CONDITION ####

Q10A_diff <- matrix(NA, ncol = 3, nrow = 7)

## Create difference in means for each condition
diff_Q10A_2 <- svyglm(Q10A_m ~ ASSIGNMENT_A_2, milconf_design)
Q10A_diff[1,1] <- diff_Q10A_2$coefficients[2]
Q10A_diff[1,2] <- confint(diff_Q10A_2)[2,1]
Q10A_diff[1,3] <- confint(diff_Q10A_2)[2,2]

diff_Q10A_3 <- svyglm(Q10A_m ~ ASSIGNMENT_A_3, milconf_design)
Q10A_diff[2,1] <- diff_Q10A_3$coefficients[2]
Q10A_diff[2,2] <- confint(diff_Q10A_3)[2,1]
Q10A_diff[2,3] <- confint(diff_Q10A_3)[2,2]

diff_Q10A_4 <- svyglm(Q10A_m ~ ASSIGNMENT_A_4, milconf_design)
Q10A_diff[3,1] <- diff_Q10A_4$coefficients[2]
Q10A_diff[3,2] <- confint(diff_Q10A_4)[2,1]
Q10A_diff[3,3] <- confint(diff_Q10A_4)[2,2]

diff_Q10A_5 <- svyglm(Q10A_m ~ ASSIGNMENT_A_5, milconf_design)
Q10A_diff[4,1] <- diff_Q10A_5$coefficients[2]
Q10A_diff[4,2] <- confint(diff_Q10A_5)[2,1]
Q10A_diff[4,3] <- confint(diff_Q10A_5)[2,2]

diff_Q10A_6 <- svyglm(Q10A_m ~ ASSIGNMENT_A_6, milconf_design)
Q10A_diff[5,1] <- diff_Q10A_6$coefficients[2]
Q10A_diff[5,2] <- confint(diff_Q10A_6)[2,1]
Q10A_diff[5,3] <- confint(diff_Q10A_6)[2,2]

diff_Q10A_7 <- svyglm(Q10A_m ~ ASSIGNMENT_A_7, milconf_design)
Q10A_diff[6,1] <- diff_Q10A_7$coefficients[2]
Q10A_diff[6,2] <- confint(diff_Q10A_7)[2,1]
Q10A_diff[6,3] <- confint(diff_Q10A_7)[2,2]

diff_Q10A_8 <- svyglm(Q10A_m ~ ASSIGNMENT_A_8, milconf_design)
Q10A_diff[7,1] <- diff_Q10A_8$coefficients[2]
Q10A_diff[7,2] <- confint(diff_Q10A_8)[2,1]
Q10A_diff[7,3] <- confint(diff_Q10A_8)[2,2]

## Plot
Q10A_plot.df <- data.frame(change = Q10A_diff[,1],
                           ci_low = Q10A_diff[,2],
                           ci_hi = Q10A_diff[,3],
                           cond = c(1:7))

Q10A_plot <- ggplot(Q10A_plot.df) + 
  geom_bar(aes(x = cond, y = change), stat = "identity", width = 0.75, fill = "gray60") +
  geom_errorbar(aes(x = cond, ymin = ci_low, ymax = ci_hi), width = 0.4)
Q10A_plot <- Q10A_plot + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1:7),
                     labels = c("2", "3", "4", "5", "6", "7", "8")) +
  ylab("Difference in Means from Control Condition") + ggtitle("Q10A by Treatment Condition")
Q10A_plot <- Q10A_plot + geom_hline(yintercept = 0, colour = "black", lty = 1)
Q10A_plot + theme_bw()

#### CONFIDENCE: Q10A BY CONDITION AND PARTY ####

Q10A_diff_p <- matrix(NA, ncol = 9, nrow = 7)

## Create difference in means, CIs for republicans
diff_Q10A_2_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_2, milconf_design)
Q10A_diff_p[1,1] <- diff_Q10A_2_r$coefficients[2]
Q10A_diff_p[1,2] <- confint(diff_Q10A_2_r)[2,1]
Q10A_diff_p[1,3] <- confint(diff_Q10A_2_r)[2,2]

diff_Q10A_3_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_3, milconf_design)
Q10A_diff_p[2,1] <- diff_Q10A_3_r$coefficients[2]
Q10A_diff_p[2,2] <- confint(diff_Q10A_3_r)[2,1]
Q10A_diff_p[2,3] <- confint(diff_Q10A_3_r)[2,2]

diff_Q10A_4_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_4, milconf_design)
Q10A_diff_p[3,1] <- diff_Q10A_4_r$coefficients[2]
Q10A_diff_p[3,2] <- confint(diff_Q10A_4_r)[2,1]
Q10A_diff_p[3,3] <- confint(diff_Q10A_4_r)[2,2]

diff_Q10A_5_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_5, milconf_design)
Q10A_diff_p[4,1] <- diff_Q10A_5_r$coefficients[2]
Q10A_diff_p[4,2] <- confint(diff_Q10A_5_r)[2,1]
Q10A_diff_p[4,3] <- confint(diff_Q10A_5_r)[2,2]

diff_Q10A_6_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_6, milconf_design)
Q10A_diff_p[5,1] <- diff_Q10A_6_r$coefficients[2]
Q10A_diff_p[5,2] <- confint(diff_Q10A_6_r)[2,1]
Q10A_diff_p[5,3] <- confint(diff_Q10A_6_r)[2,2]

diff_Q10A_7_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_7, milconf_design)
Q10A_diff_p[6,1] <- diff_Q10A_7_r$coefficients[2]
Q10A_diff_p[6,2] <- confint(diff_Q10A_7_r)[2,1]
Q10A_diff_p[6,3] <- confint(diff_Q10A_7_r)[2,2]

diff_Q10A_8_r <- svyglm(Q10A_m_rep ~ ASSIGNMENT_A_8, milconf_design)
Q10A_diff_p[7,1] <- diff_Q10A_8_r$coefficients[2]
Q10A_diff_p[7,2] <- confint(diff_Q10A_8_r)[2,1]
Q10A_diff_p[7,3] <- confint(diff_Q10A_8_r)[2,2]

## Create difference in means, CIs for independents
diff_Q10A_2_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_2, milconf_design)
Q10A_diff_p[1,4] <- diff_Q10A_2_i$coefficients[2]
Q10A_diff_p[1,5] <- confint(diff_Q10A_2_i)[2,1]
Q10A_diff_p[1,6] <- confint(diff_Q10A_2_i)[2,2]

diff_Q10A_3_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_3, milconf_design)
Q10A_diff_p[2,4] <- diff_Q10A_3_i$coefficients[2]
Q10A_diff_p[2,5] <- confint(diff_Q10A_3_i)[2,1]
Q10A_diff_p[2,6] <- confint(diff_Q10A_3_i)[2,2]

diff_Q10A_4_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_4, milconf_design)
Q10A_diff_p[3,4] <- diff_Q10A_4_i$coefficients[2]
Q10A_diff_p[3,5] <- confint(diff_Q10A_4_i)[2,1]
Q10A_diff_p[3,6] <- confint(diff_Q10A_4_i)[2,2]

diff_Q10A_5_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_5, milconf_design)
Q10A_diff_p[4,4] <- diff_Q10A_5_i$coefficients[2]
Q10A_diff_p[4,5] <- confint(diff_Q10A_5_i)[2,1]
Q10A_diff_p[4,6] <- confint(diff_Q10A_5_i)[2,2]

diff_Q10A_6_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_6, milconf_design)
Q10A_diff_p[5,4] <- diff_Q10A_6_i$coefficients[2]
Q10A_diff_p[5,5] <- confint(diff_Q10A_6_i)[2,1]
Q10A_diff_p[5,6] <- confint(diff_Q10A_6_i)[2,2]

diff_Q10A_7_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_7, milconf_design)
Q10A_diff_p[6,4] <- diff_Q10A_7_i$coefficients[2]
Q10A_diff_p[6,5] <- confint(diff_Q10A_7_i)[2,1]
Q10A_diff_p[6,6] <- confint(diff_Q10A_7_i)[2,2]

diff_Q10A_8_i <- svyglm(Q10A_m_ind ~ ASSIGNMENT_A_8, milconf_design)
Q10A_diff_p[7,4] <- diff_Q10A_8_i$coefficients[2]
Q10A_diff_p[7,5] <- confint(diff_Q10A_8_i)[2,1]
Q10A_diff_p[7,6] <- confint(diff_Q10A_8_i)[2,2]

## Create difference in means, CIs for democrats
diff_Q10A_2_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_2, milconf_design)
Q10A_diff_p[1,7] <- diff_Q10A_2_d$coefficients[2]
Q10A_diff_p[1,8] <- confint(diff_Q10A_2_d)[2,1]
Q10A_diff_p[1,9] <- confint(diff_Q10A_2_d)[2,2]

diff_Q10A_3_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_3, milconf_design)
Q10A_diff_p[2,7] <- diff_Q10A_3_d$coefficients[2]
Q10A_diff_p[2,8] <- confint(diff_Q10A_3_d)[2,1]
Q10A_diff_p[2,9] <- confint(diff_Q10A_3_d)[2,2]

diff_Q10A_4_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_4, milconf_design)
Q10A_diff_p[3,7] <- diff_Q10A_4_d$coefficients[2]
Q10A_diff_p[3,8] <- confint(diff_Q10A_4_d)[2,1]
Q10A_diff_p[3,9] <- confint(diff_Q10A_4_d)[2,2]

diff_Q10A_5_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_5, milconf_design)
Q10A_diff_p[4,7] <- diff_Q10A_5_d$coefficients[2]
Q10A_diff_p[4,8] <- confint(diff_Q10A_5_d)[2,1]
Q10A_diff_p[4,9] <- confint(diff_Q10A_5_d)[2,2]

diff_Q10A_6_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_6, milconf_design)
Q10A_diff_p[5,7] <- diff_Q10A_6_d$coefficients[2]
Q10A_diff_p[5,8] <- confint(diff_Q10A_6_d)[2,1]
Q10A_diff_p[5,9] <- confint(diff_Q10A_6_d)[2,2]

diff_Q10A_7_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_7, milconf_design)
Q10A_diff_p[6,7] <- diff_Q10A_7_d$coefficients[2]
Q10A_diff_p[6,8] <- confint(diff_Q10A_7_d)[2,1]
Q10A_diff_p[6,9] <- confint(diff_Q10A_7_d)[2,2]

diff_Q10A_8_d <- svyglm(Q10A_m_dem ~ ASSIGNMENT_A_8, milconf_design)
Q10A_diff_p[7,7] <- diff_Q10A_8_d$coefficients[2]
Q10A_diff_p[7,8] <- confint(diff_Q10A_8_d)[2,1]
Q10A_diff_p[7,9] <- confint(diff_Q10A_8_d)[2,2]

## Plot
Q10A_plot_p.df <- data.frame(rep = Q10A_diff_p[,1],
                             rep_lo = Q10A_diff_p[,2],
                             rep_hi = Q10A_diff_p[,3],
                             ind = Q10A_diff_p[,4],
                             ind_lo = Q10A_diff_p[,5],
                             ind_hi = Q10A_diff_p[,6],
                             dem = Q10A_diff_p[,7],
                             dem_lo = Q10A_diff_p[,8],
                             dem_hi = Q10A_diff_p[,9],
                             rcond = c(1,5,9,13,17,21,25),
                             icond = c(2,6,10,14,18,22,26),
                             dcond = c(3,7,11,15,19,23,27))

Q10A_plot_p <- ggplot(Q10A_plot_p.df) + 
  geom_bar(aes(x = rcond, y = rep), stat = "identity", width = 0.8, fill = "indianred1") +
  geom_errorbar(aes(x = rcond, ymin = rep_lo, ymax = rep_hi), width = 0.2) +
  geom_bar(aes(x = icond, y = ind), stat = "identity", width = 0.8, fill = "gray83") +
  geom_errorbar(aes(x = icond, ymin = ind_lo, ymax = ind_hi), width = 0.2) +
  geom_bar(aes(x = dcond, y = dem), stat = "identity", width = 0.8, fill = "steelblue2") +
  geom_errorbar(aes(x = dcond, ymin = dem_lo, ymax = dem_hi), width = 0.2)
Q10A_plot_p <- Q10A_plot_p + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6, 10, 14, 18, 22, 26),
                     labels = c("2", "3", "4", "5", "6", "7", "8")) +
  ylab("Difference in Means from Control Condition")
Q10A_plot_p <- Q10A_plot_p + geom_hline(yintercept = 0, colour = "black", lty = 1) +
  ggtitle("Q10A by Condition and Party ID")
Q10A_plot_p + theme_bw()